########################################
########################################
########################################
########################################
##Working directory should be set to location of this R file.
########################################
########################################
########################################
########################################

########################################
########################################
########################################
########################################
##Replication code for statistical analyses reported in "Legitimacy and Hegemony in Late Imperial China"
##Security Studies
##Austin Strange
##Last updated: 8 March 2024
########################################
########################################
########################################
########################################

##Load packages
require(AER)
require(data.table)
require(lmtest)
require(MASS)
require(plm)
require(stargazer)
require(xtable)
require(zoo)

########################################
########################################
########################################
########################################
##Load, merge, and clean data
########################################
########################################
########################################
######################################## 

##Load annual dataset
d1 <- read.csv("SS_AMS_ReplicationData_Annual.csv")

##Load entity-year dataset
d2 <- read.csv("SS_AMS_ReplicationData_EntityYear.csv")

##Subset for Ming and Qing separate analyses
d1m <- d1[d1$Qing == 0,]
d1q <- d1[d1$Qing == 1,]

##Subset for East Asian analysis using entity-year dataset
eas <- c("Korea", "Goryeo", "Japan", "Ryukyu")
d2ea <- d2[d2$Entity %in% eas,]

##Subset for analysis excluding Emperor Tianshun (Table 14)
d1nots <- d1[d1$Emperor != "Tianshun",]

########################################
########################################
########################################
########################################
########################################
########################################
##Replicate analysis reported in manuscript main text
########################################
########################################
########################################
########################################
########################################
########################################

########################################
########################################
##Table 1: Leader legitimacy and tribute exchanges
########################################
########################################
t1a <- lm(LExchange ~ yip.norm,
          data=d1); summary(t1a)
t1b <- lm(LExchange ~ yip.norm + factor(hcent),
          data=d1); summary(t1b)
t1c <- lm(LExchange ~ yip.norm + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB + Temp + Qing + factor(hcent),
          data=d1); summary(t1c)
t1d <- lm(LExchange ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(t1d)

t1a_rc <- coeftest(t1a,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t1b_rc <- coeftest(t1b,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t1c_rc <- coeftest(t1c,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t1d_rc <- coeftest(t1d,vcov = vcovCL,type = "HC1",cluster = ~Emperor)

stargazer(t1a_rc,t1b_rc,t1c_rc,t1d_rc,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit = c("hcent","Quart"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Period FE", "Yes", "Yes", "Yes", "Yes")),
          single.row=T)

########################################
########################################
##Table 2: Leader legitimacy, tribute exchanges, and “high-value” targets
########################################
########################################
t2a <- lm(LEA ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(t2a)
t2b <- lm(LSEA ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing +factor(hcent),
          data=d1); summary(t2b)
t2c <- lm(LIOR ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(t2c)
t2d <- lm(LNoA ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(t2d)
t2e <- lm(LCA ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(t2e)
t2f <- lm(LSW ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(t2f)

t2a_rc <- coeftest(t2a,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t2b_rc <- coeftest(t2b,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t2c_rc <- coeftest(t2c,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t2d_rc <- coeftest(t2d,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t2e_rc <- coeftest(t2e,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t2f_rc <- coeftest(t2f,vcov = vcovCL,type = "HC1",cluster = ~Emperor)

stargazer(t2a_rc,t2b_rc,t2c_rc,t2d_rc,t2e_rc,t2f_rc,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit = c("hcent"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Period FE", "Yes", "Yes", "Yes", "Yes", "Yes","Yes")),
          single.row=T)

########################################
########################################
##Table 3: Leader legitimacy and tribute (entity-year analysis)
########################################
########################################
t3a <- glm(ExB ~ yip.norm,data=d2,family = "binomial"); summary(t3a)
t3b <- glm(ExB ~ yip.norm + factor(Entity) + factor(Decade),data=d2,
           family = "binomial"); summary(t3b)
t3c <- glm(ExB ~ yip.norm + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + 
             factor(Entity) + factor(Decade),data=d2,family = "binomial"); summary(t3c)
t3d <- glm(ExB ~ yip.norm*irreg + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + 
             factor(Entity) + factor(Decade),data=d2,family = "binomial"); summary(t3d)
t3e <- glm(ExB ~ yip.norm + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + 
             factor(Entity) + factor(Decade) + L1_NLSinic,data=d2ea,family = "binomial"); summary(t3e)
t3f <- glm(ExB ~ yip.norm*irreg + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + 
             factor(Entity) + factor(Decade) + L1_NLSinic,data=d2ea,family = "binomial"); summary(t3f)

t3a_rc <- coeftest(t3a,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t3b_rc <- coeftest(t3b,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t3c_rc <- coeftest(t3c,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t3d_rc <- coeftest(t3d,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t3e_rc <- coeftest(t3e,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
t3f_rc <- coeftest(t3f,vcov = vcovCL,type = "HC1",cluster = ~Emperor)

stargazer(t3a_rc,t3b_rc,t3c_rc,t3d_rc,t3e_rc,t3f_rc,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit = c("Decade","Entity"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Decade FE", "Yes", "Yes", "Yes", "Yes", "Yes")),
          single.row=T)

########################################
########################################
########################################
########################################
########################################
########################################
##Replicate analysis reported in appendices
########################################
########################################
########################################
########################################
########################################
########################################

########################################
########################################
##Table A5: Summary Statistics for Variables Used in Main Regression Analyses
########################################
########################################
decstats <- describe(d1[,c(4,13,16,15,18:20,21:25,47,36)],skew=F,ranges = TRUE)
decstats <- setDT(decstats, keep.rownames = TRUE)[]
xtable(decstats)

########################################
########################################
#Table A6: Alternative Measures of New Leaders
########################################
########################################
a6a <- lm(LExchange ~ first2,
          data=d1); summary(a6a)
a6b <- lm(LExchange ~ first2 +factor(hcent),
          data=d1); summary(a6b)
a6c <- lm(LExchange ~ first2 + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB + Temp + Qing + factor(hcent),
          data=d1); summary(a6c)
a6d <- lm(LExchange ~ first2*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(a6d)

a6e <- lm(LExchange ~ first3,
          data=d1); summary(a6e)
a6f <- lm(LExchange ~ first3 +factor(hcent),
          data=d1); summary(a6f)
a6g <- lm(LExchange ~ first3 + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB + Temp + Qing + factor(hcent),
          data=d1); summary(a6g)
a6h <- lm(LExchange ~ first3*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(a6h)

a6a_rc <- coeftest(a6a,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a6b_rc <- coeftest(a6b,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a6c_rc <- coeftest(a6c,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a6d_rc <- coeftest(a6d,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a6e_rc <- coeftest(a6e,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a6f_rc <- coeftest(a6f,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a6g_rc <- coeftest(a6g,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a6h_rc <- coeftest(a6h,vcov = vcovCL,type = "HC1",cluster = ~Emperor)

stargazer(a6a_rc,a6b_rc,a6c_rc,a6d_rc,a6e_rc,a6f_rc,a6g_rc,a6h_rc,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit = c("hcent"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Period FE", "Yes", "Yes", "Yes", "Yes", "Yes","Yes", "Yes", "Yes")),
          single.row=T)

########################################
########################################
#Table A7: Alternative Measures of Irregular Entry
########################################
########################################
a7a <- lm(LExchange ~ yip.norm*Coup+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(a7a)
a7b <- lm(LExchange ~ yip.norm*irreg2+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(a7b)

a7a_rc <- coeftest(a7a,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a7b_rc <- coeftest(a7b,vcov = vcovCL,type = "HC1",cluster = ~Emperor)

stargazer(a7a_rc,a7b_rc,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit = c("hcent"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Period FE", "Yes", "Yes")),
          single.row=T)

########################################
########################################
#Table A8: Leader Legitimacy and Tribute Exchanges (Negative Binomial)
########################################
########################################
a8a <- glm.nb(Exchange ~ yip.norm,
          data=d1); summary(a8a)
a8b <- glm.nb(Exchange ~ yip.norm + factor(hcent),
          data=d1); summary(a8b)
a8c <- glm.nb(Exchange ~ yip.norm + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB + Temp + Qing + factor(hcent),
          data=d1); summary(a8c)
a8d <- glm.nb(Exchange ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(a8d)

a8a_rc <- coeftest(a8a,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a8b_rc <- coeftest(a8b,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a8c_rc <- coeftest(a8c,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a8d_rc <- coeftest(a8d,vcov = vcovCL,type = "HC1",cluster = ~Emperor)

stargazer(a8a_rc,a8b_rc,a8c_rc,a8d_rc,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit = c("hcent"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Period FE", "Yes", "Yes", "Yes", "Yes")),
          single.row=T)

########################################
########################################
#Table A9: Leader Legitimacy and Tribute Exchanges (Raw annual counts)
########################################
########################################
a9aCount <- lm(Exchange ~ first.third,
               data=d1); summary(a9aCount)
a9bCount <- lm(Exchange ~ first.third + factor(hcent),
               data=d1); summary(a9bCount)
a9cCount <- lm(Exchange ~ first.third +  L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB + Temp + Qing + factor(hcent),
               data=d1); summary(a9cCount)
a9dCount <- lm(Exchange ~ first.third*irreg+  L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
               data=d1); summary(a9dCount)

a9a_rcCount <- coeftest(a9aCount,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a9b_rcCount <- coeftest(a9bCount,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a9c_rcCount <- coeftest(a9cCount,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a9d_rcCount <- coeftest(a9dCount,vcov = vcovCL,type = "HC1",cluster = ~Emperor)

stargazer(a9a_rcCount,a9b_rcCount,a9c_rcCount,a9d_rcCount,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit = c("hcent"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Period FE", "Yes", "Yes", "Yes", "Yes")),
          single.row=T)

########################################
########################################
#Table A10: Leader Legitimacy and Tribute Exchanges (Table 1 w/ all covariates)
########################################
########################################
a10a <- lm(LExchange ~ yip.norm,
          data=d1); summary(a10a)
a10b <- lm(LExchange ~ yip.norm + factor(hcent),
          data=d1); summary(a10b)
a10c <- lm(LExchange ~ yip.norm + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB + Temp + Qing + factor(hcent),
          data=d1); summary(a10c)
a10d <- lm(LExchange ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(a10d)

a10a_rc <- coeftest(a10a,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a10b_rc <- coeftest(a10b,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a10c_rc <- coeftest(a10c,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a10d_rc <- coeftest(a10d,vcov = vcovCL,type = "HC1",cluster = ~Emperor)

stargazer(a10a_rc,a10b_rc,a10c_rc,a10d_rc,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit = c("hcent","Quart"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Period FE", "Yes", "Yes", "Yes", "Yes")),
          single.row=T)

########################################
########################################
#Table A11: Leader Legitimacy, Tribute Exchanges, and “High-Value” Targets (Table 2 w/ all Covariates)
########################################
########################################
a11a <- lm(LEA ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(a11a)
a11b <- lm(LSEA ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing +factor(hcent),
          data=d1); summary(a11b)
a11c <- lm(LIOR ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(a11c)
a11d <- lm(LNoA ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(a11d)
a11e <- lm(LCA ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(a11e)
a11f <- lm(LSW ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
          data=d1); summary(a11f)

a11a_rc <- coeftest(a11a,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a11b_rc <- coeftest(a11b,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a11c_rc <- coeftest(a11c,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a11d_rc <- coeftest(a11d,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a11e_rc <- coeftest(a11e,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a11f_rc <- coeftest(a11f,vcov = vcovCL,type = "HC1",cluster = ~Emperor)

stargazer(a11a_rc,a11b_rc,a11c_rc,a11d_rc,a11e_rc,a11f_rc,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit = c("hcent"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Period FE", "Yes", "Yes", "Yes", "Yes", "Yes","Yes")),
          single.row=T)

########################################
########################################
#Table A12: Leader Legitimacy and Tribute, Entity-Year Analysis (Table 3 w/ all covariates)
########################################
########################################
a12a <- glm(ExB ~ yip.norm,data=d2,family = "binomial"); summary(a12a)
a12b <- glm(ExB ~ yip.norm + factor(Entity) + factor(Decade),data=d2,
           family = "binomial"); summary(a12b)
a12c <- glm(ExB ~ yip.norm + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + 
             factor(Entity) + factor(Decade),data=d2,family = "binomial"); summary(a12c)
a12d <- glm(ExB ~ yip.norm*irreg + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + 
             factor(Entity) + factor(Decade),data=d2,family = "binomial"); summary(a12d)
a12e <- glm(ExB ~ yip.norm + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + 
             factor(Entity) + factor(Decade) + L1_NLSinic,data=d2ea,family = "binomial"); summary(a12e)
a12f <- glm(ExB ~ yip.norm*irreg + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + 
             factor(Entity) + factor(Decade) + L1_NLSinic,data=d2ea,family = "binomial"); summary(a12f)

a12a_rc <- coeftest(a12a,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a12b_rc <- coeftest(a12b,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a12c_rc <- coeftest(a12c,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a12d_rc <- coeftest(a12d,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a12e_rc <- coeftest(a12e,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a12f_rc <- coeftest(a12f,vcov = vcovCL,type = "HC1",cluster = ~Emperor)

stargazer(a12a_rc,a12b_rc,a12c_rc,a12d_rc,a12e_rc,a12f_rc,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit = c("Decade","Entity"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Decade FE", "Yes", "Yes", "Yes", "Yes", "Yes")),
          single.row=T)

########################################
########################################
#Table A13: Leader Legitimacy and Tribute (Ming and Qing)
########################################
########################################
a13c <- lm(LExchange ~ yip.norm +  L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB + Temp + factor(hcent),
              data=d1m); summary(a13c)
a13d <- lm(LExchange ~ yip.norm*irreg+  L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + factor(hcent),
              data=d1m); summary(a13d)
a13g <- lm(LExchange ~ yip.norm +  L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB + Temp + factor(hcent),
          data=d1q); summary(a13g)
a13h <- lm(LExchange ~ yip.norm*irreg+ L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + factor(hcent),
          data=d1q); summary(a13h)

a13c_rc <- coeftest(a13c,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a13d_rc <- coeftest(a13d,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a13g_rc <- coeftest(a13g,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a13h_rc <- coeftest(a13h,vcov = vcovCL,type = "HC1",cluster = ~Emperor)

stargazer(a13c_rc,a13d_rc,a13g_rc,a13h_rc,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit = c("hcent"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Period FE", "Yes", "Yes", "Yes", "Yes")),
          single.row=T)

########################################
########################################
#Table A14: Main Analysis (Table 1) without Emperor Tianshun
########################################
########################################
a14a <- lm(LExchange ~ yip.norm,
             data=d1nots); summary(a14a)
a14b <- lm(LExchange ~ yip.norm + factor(hcent),
             data=d1nots); summary(a14b)
a14c <- lm(LExchange ~ yip.norm + L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB + Temp + Qing + factor(hcent),
             data=d1nots); summary(a14c)
a14d <- lm(LExchange ~ yip.norm*irreg+ L1_wokouB + L1_uprisingB + L1_JapanB + L1_maritimeB + L1_impmilB  + Temp + Qing + factor(hcent),
             data=d1nots); summary(a14d)

a14a_rc <- coeftest(a14a,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a14b_rc <- coeftest(a14b,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a14c_rc <- coeftest(a14c,vcov = vcovCL,type = "HC1",cluster = ~Emperor)
a14d_rc <- coeftest(a14d,vcov = vcovCL,type = "HC1",cluster = ~Emperor)

stargazer(a14a_rc,a14b_rc,a14c_rc,a14d_rc,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit = c("hcent"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Period FE", "Yes", "Yes", "Yes", "Yes")),
          single.row=T)